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Abstract 

This paper concerns the reconstruction of a diffusion coefficient in an ehiptic 
equation from knowledge of several power densities. The power density is the prod- 
uct of the diffusion coefficient with the square of the modulus of the gradient of the 
elliptic solution. The derivation of such internal functionals comes from perturbing 
the medium of interest by acoustic (plane) waves, which results in small changes 
in the diffusion coefficient. After appropriate asymptotic expansions and (Fourier) 
transformation, this allow us to construct the power density of the equation point- 
wise inside the domain. Such a setting finds applications in ultrasound modulated 
electrical impedance tomography and ultrasound modulated optical tomography. 

We show that the diffusion coefficient can be uniquely and stably reconstructed 
from knowledge of a sufficient large number of power densities. Explicit expressions 
for the reconstruction of the diffusion coefficient are also provided. Such results 
hold for a large class of boundary conditions for the elliptic equation in the two- 
dimensional setting. In three dimensions, the results are proved for a more restrictive 
class of boundary conditions constructed by means of complex geometrical optics 
solutions. 

1 Introduction 

Optical tomography (OT) and electrical impedance tomography (EIT) are medical imag- 
ing techniques that utilize the large contrast between the optical and electrical response 
of certain unhealthy tissues and that of healthy tissues. However, because OT and EIT 
are modeled by diffusion operators that are highly smoothing, the reconstructions are 
typically very low resolution [3l [5l [27] . Several recent imaging techniques aim to combine 
the high contrast in OT and EIT with a high resolution modality, for instance based on 
magnetic resonances (MREIT) [191 EH [22] or, for what is of interest for us here, ultra- 
sound. There are several ways of combining EIT and OT with ultrasound. One way is 
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to use the photo-acoustic or the electro-acoustic effects; see, e.g., [HI [121 [221 [2B]- In this 
paper, we consider another way based on ultrasound modulation. 

In the setting of EIT, ultrasound modulated electrical impedance tomography (UMEIT) 
is also called acousto-electric tomography; see [21 [6l [121 [El [HI [29]. In the setting of opti- 
cal tomography, ultrasound modulated tomography (UMOT) is also called acousto-optic 
tomography. A model has been derived in [S] in the so-called incoherent regime of wave 
propagation, while the coherent regime is addressed in, e.g., |^ \T7\ [28]. 

The inverse problem of UMEIT may be described as follows. The electric potential u 
solves a diffusion equation (see ([T|) below) with unknown conductivity a. The ultrasound 
modulation whose modeling is recalled in section ([2]) provides access to the power density 
H{x) = a{x)\'Vu\'^{x) for all x inside the domain of interest. Reconstruction of a from 
this one internal functional H{x) is typically not feasible unless additional constraints 
are imposed; see e.g., [6]. In this paper, we assume that several such power densities 
corresponding to, e.g., different boundary conditions in ([1]) are known. By polarization 
or direct measurements as shown in section [21 we can have access to Hij{x) = a{x)'Vui ■ 
S/Uj{x) for Ui and Uj solutions to the same elliptic equation with different boundary 
conditions. The case of the spatial dimension n = 2 was treated in [T2|. This paper aims 
to generalize these results to the case n = 3 (which also generalize to the case n > 4, see 
|20] . although we shall not present the details here) and to show that UMEIT is a stable 
inverse problem, in the sense that in an appropriate norm, errors on the measurement of 
the power densities Hij result in errors on the reconstruction of a that are of the same 
order. This should be contrasted with the case of EIT and OT, where the error on cr 
is roughly proportional to the logarithm of (and hence much larger than) the error on 
the available (Cauchy) data. Moreover, our derivation is constructive in the sense that 
it produces an explicit procedure to reconstruct a from knowledge of a sufficiently large 
number of internal measurements of the form Hy. 

In dimension n = 2, one set of measurements for 1 < i, j < 2 is sufficient to uniquely 
characterize a. In dimension n = 3, sets of measurements 1 < i, j < 3 are sufficient locally 
but not necessarily globally. We show the existence of sets of measurements 1 < j < 4, 
which uniquely (and stably) characterize a. This construction is based on showing that 
the determinant of a matrix whose columns are the gradients of three of the four solutions 
to the elliptic equation is always positive at any point x of the domain. The three chosen 
solutions depend on the point x. Such constructions are based on the use of complex 
geometric optics solutions, which are a convenient tool to show that solutions to elliptic 
equations satisfy desired qualitative properties [6l [3 [HI [26] . 

The rest of the paper is structured as follows. The derivation of the power density 
measurements from ultrasound modulations of media of interest is recalled in section [21 
The main results obtained in this paper are presented in section [21 The proofs of the main 
theorems are given in the following two sections. Section [D provides the proof of unique- 
ness, stability, and reconstruction procedures on a domain Q on which the determinant 
of a matrix composed of n gradients of solutions is shown to be bounded below by a posi- 
tive constant. How such results can be patched together to provide global reconstruction 
results is given in section [3 That determinants can be shown to have a prescribed sign 
is based on the use of complex geometric optics (CGO) solutions that are presented in 
section [51 
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2 Derivation of internal functionals 



We consider the following elliptic equation 



V ■ cr(x)Vu = in X, 



u = g on dX. 



(1) 



Here, a is the unknown diffusion coefficient, which we assume in this paper is a real- 
valued, scalar, function defined on a simply connected domain X C M" for n = 2 or 
n = 3. We assume that a is bounded above and below by positive constants so that the 
above equation admits a unique solution. We also assume that a is sufficiently smooth 
so that the solution to the above equation is continuously differentiable on X, the closure 
of X [15]. We denote by dX the boundary of X and by g{x) the imposed (sufficiently 
smooth) Dirichlet boundary conditions. 

The coefficient a{x) may model the electrical conductivity in the setting of electrical 
impedance tomography (EIT) or a diffusion coefficient of particles (photons) in the setting 
of optical tomography (OT). Both EIT and OT are modalities with high contrast, in the 
sense that (t{x) takes different values in different tissues and allows one to discriminate 
between healthy and non-healthy tissues. In OT, high contrasts are mostly observed in 
the absorption coefficient, which is not modeled here |H]. 

A methodology to couple high contrast with high resolution consists of perturbing the 
diffusion coefficient acoustically. Let an acoustic signal propagate through the domain. 
In this presentation, we assume that the sound speed is constant and that the acoustic 
signal is a plane wave of the form pcos{k ■ x + (f) where p is the amplitude of the acoustic 
signal, k its wave-number and (f an additional phase. The acoustic signal modifies the 
properties of the diffusion equation. We assume that such an effect is small and that the 
coefficient in ([1]) is modified as 



where we have defined c = c{x) = cos(/c ■ x + ip) and where e = pT is the product of 
the acoustic amplitude p G M and a measure F > of the coupling between the acoustic 
signal and the modulations of the constitutive parameter in ([1]), see [18] and references 
therein. We assume that e ^ 1 so that the influence of the acoustic signal on admits 
an asymptotic expansion that we truncated at the second order as displayed in ([2]). The 
size of the terms in the expansion are physically characterized by ( and depend on the 
specific application. 

Let u and v be solutions of ([1]) with fixed boundary conditions g and h, respectively. 
When the acoustic field is turned on, the coefficients are modified as described in ([2]) 
and we denote by and the corresponding solutions. Note that M_e is the solution 
obtained by changing the sign of p or equivalently by replacing ip hj ip + it. 

By the standard continuity of the solution to ([1]) with respect to changes in the coef- 
ficients and regular perturbation arguments, we find that = uq + eui + 0(6^). Let us 
multiply the equation for by and the equation for v^^ by u^, subtract the resulting 
equalities, and use standard integrations by parts. We obtain that 




(2) 



I ((Te - a_e)Vn£ ■ Vv_^dx 

X 




dUe 



V-sds{x). 



(3) 



du 



dv 
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Here, ds{x) is the standard surface measure on dX. We assume that a^dyUi; and a^d^Vi; 
are measured on dX, at least on the support of Vs = h and Us = g, respectively, for 
all values e of interest. Note that the above equation holds if the Dirichlet boundary 
conditions are replaced by Neumann boundary conditions. Let us define 

1 f Ov du 

J, :=-/ a^e^^u,-a,—^v^eds{x)=e.h+e^J2 + 0{e'^). (4) 

2 J ax ov du 

We assume that the real valued functions = Jm{k, ^) are known (measured functions). 
Notice that such knowledge is based on the physical boundary measurement of the Cauchy 
data of the form (m^, aed^Ue) and (wg, a^d^v^) on dX. 

Equating like powers of e, we find that at leading order 

/[C.(x)V„,„.V„„W]cos(«:.x + ^)dx = J,(«:.^). (5) 

Jx 

This may be acquired for all G M" and ip = 0, ^, and hence provides the Fourier 
transform of 

H[uo,vo]{x) = Co-(x)Vno ■ Vvo{x). (6) 

Note that when = Me, then we find from the expression in ([3]) that J2 = in (jlj) so 
that the expression for Ji may be obtained from available measurements in (j4]) with an 
accuracy of order O(e^). Note also that 

H[uo,Vo]{x) = ^{H[uo + Vo,Uo + Vo] - H[uo - Vo,Uo - Vq]) 

by polarization. In other words, the limiting measurements (for small e) in (E]) may also 
be obtained by considering expressions of the form ([3]) with Us = Ve- 
in the setting of optical tomography, the coefficient in ([2]) takes the form 



where is the diffusion coefficient, q is the light speed, and n is spatial dimension. When 
the pressure field is turned on, the location of the scatterers is modified by compression 
and dilation. Since the diffusion coefficient is inversely proportional to the scattering 
coefficient, we find that 

-^ = ^{l+ec{x)). 
as{x) a[x) 

Moreover, the pressure field changes the index of refraction (the speed) of light as follows 

Ce{x) = c(x)(l + 7£:c(x)), 
where 7 is a constant (roughly equal to | for water). This shows that 

C = -(l + (rf-l)7). (7) 

In the setting of electrical impedance tomography, we simply assume that ( models the 
coupling between the acoustic signal and the change in the electrical conductivity of the 
underlying material. The value of ( thus depends on the application. 
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The objective of ultrasound modulated optical tomography (UMOT) and ultrasound 
modulated electrical impedance tomography (UMEIT) is to reconstruct the coefficient 
(t{x) from measurements of the form ([6]), i.e., since we assume that ( is known, from 
measurements of the form H{x) = a{x)'Vu{x) ■ Vv{x), where u and v are two solutions of 
the unperturbed equation with (possibly) different boundary conditions. 

In the setting where Vq = Uq, measurements are of the form Hqo{x) = cr(x)|VMop. 
Plugging the latter expression into the elliptic equation yields the nonlinear equation 

V-^^Vuo = inX, uo = g ondX. (8) 

We thus observe that the reconstruction of a may be recast as solving the above nonlinear 
partial differential equation. It turns out that the above equation is not directly amenable 
to analysis. However, the above expression makes clear the connection between inverse 
problems with internal information and nonlinear partial differential equations. 

The methodology to obtain ([6]) follows the presentation in [8] and is very similar in 
spirit to the derivation obtained in [18] . An alternative method based on physical focusing 
in time of acoustic signals has been presented in [2]. 

Assuming the validity of the above derivation, the mathematical problem of interest 
in this paper is as follows. We want to reconstruct cr(x) from knowledge of the interior 
functionals 

Hij(x) = (T{x)'Vui{x) ■ Vuj{x), i < i, j < m, (9) 
where Uj is the solution to the equation 

' (10) 
Ui = Qi oX, 1 < i < m, 

for appropriate choices of the boundary conditions Qi on dX. Obviously, we would like m 
to be as small as possible and ideally equal to 1. Our results propose explicit reconstruc- 
tions for m = n in dimension n = 2 and m = n or m = n + l depending on additional 
parameters in dimension n = 3. The case of the dimension n > 4 can be handled in 
a similar fashion although we shall not present the mathematically more involved and 
practically less interesting details in this paper. 



3 Statement of the main results 

The main strategy we use to reconstruct a from knowledge of H = {Hij}ij is ffist to 
write equations for the quantities Si := y/aVui that are independent of a and then to 
show that a is uniquely determined when Si is known. This was implemented in the two 
dimensional setting in [12] and the main result of this paper is the extension to the three 
dimensional case. The extension to arbitrary dimensions is also feasible although we shall 
restrict ourselves to the cases n = 2, 3 for concreteness. 
We define: 

^. := y^Vn,, l<i<m, F := V(log ^a) = loga. (11) 
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Using the equations ffTUi) . the definitions ffTT]) and the fact that a ^ Si is a gradient, we 
obtain: 



V • iV^Sj) = ^ V ■ Sj + F ■ Sj = 0, (12) 



n = 2 : 



n 



^ [V,S',]-[F,5',] = 0, (13) 
3 : V X (^•S'j^ =0 ^ V X Sj - F X Sj = 0, (14) 



where we have defined for n = 2 the product [A,B] := A^By — AyB^ and [V,A] := 
dxAy — dyAx for smooth vector fields A and 5, while for n = 3, x is the standard cross 
product. 

We now wish to eliminate F from such equations and get a closed form equation for 
the vectors Si with sources that only involve the known matrix H . Such an elimination 
requires that we find n vectors Si that form a basis of M" for = 2, 3. In dimension 
n = 2, it is not difficult to find such a basis for all x G X. In dimension n > 3, we are 
able to construct such bases on subset of M" that cover X using triplets of vectors that 
depend on the subset. We thus consider the above equations for Si and F over an open 
subset f2 C X, over which we make the assumption that 

infdet(5i(x),...,5„(x)) = co>0, (15) 

where the n vectors Sj for 1 < j < n are chosen among the m vectors considered in (|9]). 
Since the determinants are central in our derivations, we define 

(i(x) := det(S'i(x), S'2(x)), n = 2 and D{x) := det{Si{x), S2{x), S^i^x)), n = 3, (16) 

for all X en. Note that det(^i(x), . . . , Snix)) = {detH)^x) > Cq on fi. 

Remark 3.1. Since H is invertible and is a symmetric non-negative matrix, it satisfies 
the following inequalities 



|i7-^||-"<deti/<n^.^< ll^l 



i=l 

where the second inequality is a general property of gramian matrices (i.e. matrices of 
dotproducts of a given family of vectors). Indeed, for n = 2, we have detH = H11H22 — 
^ H11H22, and for n = 3, we may write 

det H = (det(5i, 52, S-,))' = \\S, ■ {S2 x Ss)f < \\S,f\\S2 x S^f < \\S,f\\S2f\\Ssf, 

after successively using the Cauchy-Schwarz inequality and a property of the cross-product. 
Then, since H = S'^S G C^{Q) (by regularity of the solutions and assuming a continuous), 
the assumption ( fTSl) is equivalent to H being invertible over Vt. In this case is 
uniformly bounded over Q, that is ||iJ"^||oo < +00, and so 

0<\\H-'U<cl. (17) 
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In two dimensions, the constraint fll5p is satisfied over the whole domain Q = X ioi 
a large class of boundary conditions Qi and we then choose m = n = 2 as shown in 
[H Theorem 4]. For instance, we may choose gj{xi,X2) = Xj on dX and we are then 
guaranteed that (xi,X2) i— )■ (mi,M2) is a diffeomorphism from X to its image so that the 
determinant is clearly signed, and by continuity (assuming that a is sufficiently smooth 
so that gradients are continuous functions, which we also assume for the rest of the paper 
[15]), bounded from below by a positive constant. 

In three dimensions, there is no known guarantee that the determinant of three so- 
lutions with prescribed boundary conditions remains positive throughout the domain X 
independent of the conductivity a. For boundary conditions of the form gj{x) = Xj on dX 
for 1 < j < 3, it is known that the determinant can change signs for some conductivities 

m- 

Yet, (fT5|) is an important constraint in the derivation of our results. By means of 
specific complex geometric optics solutions, we are able to construct boundary condi- 
tions Qi such that f llSp is valid on subsets of the domain X. More precisely, we shall 
construct gfj for 1 < i < 4 so that det(S'i(x), S'2(x), S'3(x)) > Cq > on U^^^^Q2k and 
det(S'i(a:), 5'2(x), S'4(x)) > Cq > on U^^^ifl2k-i with all domains Qk simply connected, 
open and such that X C Ul'^f^Qk- In other words, we can choose m = 4 boundary con- 
ditions such that X is decomposed into a superposition of overlapping simply connected 
subsets where the determinant of three out of the four vectors is positive. Generically, 
we denote by Q any of the subdomains Q^- We shall see that a can be reconstructed on 
all of Q provided that a and all Si are known at one point Xq G Q. After patching local 
reconstructions together, this provides a (unique) reconstruction in the whole domain X. 

On each of the subdomains Q = Qk, the available information is Hij{x) = Si{x) ■ Sj{x) 
for 1 < i < j < n, abbreviated by a function H : Q ^ Sn(M.) (S'„(R) denotes the n x n 
real- valued symmetric matrices). Although our goal is ultimately to recover only the 
conductivity a, our approach requires the reconstruction of the vector fields Si, 1 < i < m 
and F. Since the data H{x) give us S{x) (the matrix with columns given by the vectors 
Si) up to an ^(^^(ffi)- valued function R{x) (i.e., R{x) is a rotation matrix), the unknowns 
are now the functions R and F, which are of dimension and n, respectively. As 

will be given in more detail in section HJ the local reconstruction of R and F proceeds 
as follows: one first derives divergence and curl equations for the Rj^s, column vectors of 
the function i? : f2 — )■ S'0„(M), with right-hand-sides involving F, the Rj^s and the data. 
Using structural properties of ^(^^(IR) satisfied at every x E fl, one is then able to derive 
an equation for F of the form 

i^ = M^Vlogdeti7+ J] ((\/,, + V^,,)-^.)^. ) , n = 2,3, (18) 

where the vector fields Vij depend only on the data. Plugging this equation back into 
the divergence/curl system closes the system for the i?j's. From this closed sytem, one is 
then able to derive gradient-type equations for the i?j's, where the right-hand sides either 
depend only on the data (case = 2), or also depend polynomially on the Rj's (case 
n = 3). In either case, these gradient equations can be integrated along segments of the 
form [xo,x] for x,Xo G Q, parameterized by the curve 

7xo,x : [0, 1] 9 t ^ 7xo,x(t) = (1 - t)xo + tx E n, (19) 
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in order to reconstruct locally the i?j's at x from knowledge of their value at xq. Once the 
i?j's (and therefore the function R) are locally reconstructed around xq, one can recon- 
struct (t(x) around xq from the knowledge of o"(xo) and integrating (ITSl) along segments 
[xo,x]. 

As we mentioned earlier, we can set = X when n = 2 by choosing illuminations 
g = {giiQi) such that (ITS]) is satisfied throughout X. Then the above reconstruction 
procedure yields unique and stable reconstructions, as described in the following theorem. 

Theorem 3.2 (2D global uniqueness and stability). Let {H,H) be two data sets with 
coefficients in W^'°°{X) corresponding to the same illumination g = {gi,g2) and with 
determinants satisfying the condition 

inf (rf, d)>co> 0, with d^{x) = det H{x), (P{x) = det H{x). (20) 

Let a and a be the corresponding conductivities and assume that cr{xo) = a{xo) for some 
Xo G X. Then, we have the following stability estimate 

II logo- - log5-||vyi.oo(jf) < C\\H - /f||,4/i,oo(x)- (21) 

This result also ensures uniqueness, since ( I2T1) implies H = H =^ a = a . 

Note that theorem 13.21 does not require the prescription of the S'j's at the point Xq 
because the data H together with the maximum principle allow us to access their values at 
appropriate boundary points, as will be seen in section 14.2. paragraph "reconstruction 
procedure" . 

In the three dimensional mentioned earlier, we do not know of any illumination 

g that guarantees the constraint f lT5|) throughout X. We then work with m > 3 solutions 
of ( fTOi) and assume that there exists an open covering O = {^k}i<k<N {X C U^^^flk), a 
constant Cq > and a function t : [1, N] 3 i T{i) = {T{i)i,T{i)2,T{i)s) G [l,m]^, such 
that 

inf det{Sr(i)^{x), 5*^(^)2 (x), 5*^(^)3 (x)) > cq, 1 < i < iV. (22) 

The next lemma, whose proof is given in the beginning of section [5l gives an example 
of such a setting based on the use of complex geometric optics solutions, if one assumes 
some regularity on a. 

Lemma 3.3. Let n = 3 and a G H^~^^{X) for some e > 0, be bounded from below by a 
positive constant. Then there exists a non-empty open set G C {H^i^dX))^ of quadruples 
of illuminations such that for any g = {gi, g2, 93, 94,) G G, there exists an open cover of X 
of the form {Q2i~i, ^2i}i<i<N o,nd a constant cq > such that 

inf det (5*1, 5*2, Q 5*4) > cq and inf det(S'i, 6*2, e^Ss) > cq, 1 < « < iV, (23) 

for ei and q equal to ±1. 

Using m > 3 solutions that satisfy assumption f l22p . the local reconstruction approach 
can be applied to reconstruct a on each of the fij's, and some additional work is then 
required to patch these reconstructions together and show that the global reconstruction 
scheme ensures unique and stable reconstructions of a over X. Global reconstruction 
procedures are thus presented in section [5] and summarized in the following: 
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Theorem 3.4 (3D global uniqueness and stability). Let X (Z M.^ be an open convex 
bounded domain, and let two sets of m > 3 solutions of (fTOj) generate measurements 
{H,H) whose components belong to W^'°°{X). Assume that one can define a couple 
{O, t) such that (122!) is satisfied for both sets of solutions S and S. Let also xq G Vti^ C X 
anda{xo), cr(xo), {Sr{io)i{xo), Sr(io)X^o)}i<i<3 be given. Let a and a be the conductivities 
corresponding to the measurements H and H , respectively. Then we have the following 
stability estimate: 

II logcr - \oga\\w^,^{x) < C{eo + \\H - H\\wi-^c^(x)) , (24) 
where eo is the error at the initial point Xo 

3 

eo = I log (To -log ao I + ^\\Sr{io)Xxo) - '^r{io),(a;o)||- 

1=1 

The uniqueness and stability estimate results are constructive. As we shall see, the 
reconstruction depends on the choice of curves along which the sets of ordinary differential 
equations coming from the gradient equations are solved. In the presence of noise-free 
data, the reconstruction is unique as we just presented. In the presence of noise, the 
available data may then no longer be in the range of the measurement operator mapping 
the unknown coefficient a to the measurements H{x). In this case, the reconstruction 
may depend on the choice of the curve along which integrations are performed. The 
compatibility conditions that the data need to verify in order for the reconstruction to be 
independent of the choice of path (Poincare-type lemma) is straightforward in dimension 
n = 2 but seems to be more challenging in dimension n = 3. 

Both theorems provide us with Lipschitz constants of stability in dimensions n = 2 
and n = 3. These estimates show that reconstructions of the conductivity from knowledge 
of a sufficient number of power densities is a well-posed problem, unlike what is observed 
in the Calderon problem, which consists of reconstructing the diffusion coefficient from 
knowledge of the Dirichlet to Neumann map ^7\. Thus UMEIT and UMOT allow us to 
combine the high contrast of optical and electrical properties of domains with the high 
resolution capabilities of ultrasound. 

The rest of the paper is structured as follows. In section |4] we derive the formulas that 
are necessary for local (global when n = 2) reconstructions and stability results. Section 
E] explains the global reconstruction scheme in the case n = 3 and includes the proof of 
theorem 13.41 

4 Derivation of local reconstruction formulas 

In this section, we prove Theorem 13.21 and Theorem 13.41 in the setting where X is replaced 
by an open bounded subset f2 C X C such that fllSp holds. The proof of Theorem 13.41 
in the case of arbitrary X will be concluded in section O 

Knowledge of the matrix H{x) together with the determinant condition f|T5|) allow us 
to reconstruct the matrix-valued function S{x) := [5'i(x)| . . . |5'„(x)] up to an S'0„(M)- 
valued function R{x) = [-Ri(a;)| . . . \Rn{x)]. This fact can be seen for instance by noticing 
that the orientation-preserving Gram-Schmidt procedure that creates the orthonormal 
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i?j's from the S'j's with det-RdetS' > only involves coefficients that depend on the 
inner products Si ■ Sj = Hij, which are indeed known. Note that SH~^{x) is also a 
rotation- valued function on Q. 

More generally, we take T{x) = {tijix)}^^^ a matrix-valued function on Q that 
satisfies the properties 

T^{x)T{x) = H~\x) and < detT(x) = (detif(x))"5, x E n. (25) 

Two examples are the symmetric T = and the lower-triangular T in (!59l) below 

obtained by the Gram Schmidt procedure. Denote T the matrix constructed from H in 
the same manner as T. We impose the existence of a constant Ct > such that 

\\T — r||vKi.°°(Q) < Ct\\H — if||vFi.°°(f7), (26) 

here Ct only depends on H and H and is completely determined by the way we construct 
T and T. The following simple lemma, whose proof is given in section I4.2.3[ justifies the 
above stability result for the Gram-Schmidt procedure: 

Lemma 4.1. Let H,H E W^'°°{Q) that are invertible over Q and set 

mo = max "^||oo, ^\\oo^ , Mq = max (^\\H ^\\oo\\H\\wi,^(^q^^\\H ^\\oo\\H\\wi'°°{n)j ■ 

If (T, T) are built from {H, H) using the Gram-Schmidt procedure, then the stability esti- 
mate holds 

\\T - T||H/i.oo(n) < ml^''Cn{Mo)\\H - H\\wu..^n), (27) 
where C„(^) is a real polynomial function in C, that only depends on the dimension n. 

We also define T^^ = {t''^}i<ij<n and the vector fields 

:= V(t,fc)t'^ i.e. := 9Kt.fc)t'^ l<t,jj< n. (28) 

Here and below, repeated indices are summed over. Let V he a vector-valued matrix 
constructed from if as is constructed from H. A simple calculation yields 

\\Vij - Vij\\L^(^n) < Cv\\H - H\\wi,c^^Q), l<i,j <n. (29) 

Here, Cy = (||i/-i^'||ii||oo + ||^-i^'||^||oo) Ct + ||^~ioo, with Ct defined in ([26]). 

From 5* and T, let us now build the matrix R by defining 

R{x) := S{x)T{xf, i.e. R^{x) =tij{x)Sj{x), 1 < i < n, x E n. (30) 

From conditions f l25|) and the fact that S{x)'^S{x) = H{x) for every x E Q, the matrix R 
(pOj) satisfies R^{x)R{x) = I„ for all x G fi, as well as det R{x) = 1, thus R{x) E 50„(M) 
for all X E Q. Moreover, plugging relation f l30|) into equation f fT2|) and using the vector 
calculus identity V ■ (fV) = V/ ■ V -\- /V ■ V, we obtain 

V ■ = V ■ {t.jSj) = {Vtij) ■ Sj + tijV ■ Sj = iVt,,)f'' ■ Rk - tijF ■ S,. 
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Thus the -Rj's satisfy the following divergence equation: 

V ■ i?. = V,k -Rk-F-R,, l<i<n. (31) 

In a similar manner, one can derive the following curl-type equations for the -Rj's: 

n = 2 : [V, R^] = [Vik, Rk] + [F, Ri], i = 1, 2, (32) 
n = 3 : V X Ri = V^kX Rk + F X Ri, i = 1,2,3. (33) 

We now show that the redundancies in the systems (|3T])- (!32|) when n = 2 and ( 13T|) - (133|) 
when n = 3 allow us to derive formula (fT8|) for F in terms of the i?j's and the data. This 
formula is successively proved for n = 2 and n = 3 using vector calculus identities in each 
dimension, respectively. 

Remark 4.2 (Geometric motivation of the next proofs.). The generalization to n > 4 
may be found in I20f . where the natural tool is differential geometry, in particular connec- 
tions, exterior derivatives and the Hodge transform. Once condition (IT^ is assumed, an 
expression of equation (ITS]) in the basis {Si, . . . , Sn) may be proved in general dimension 
by studying the properties of the dual basis of [Si, . . . , Sn) with respect to the Euclidean 
metric. Once f|T8|) is derived, the divergence and curl equations that we have for the Si 's 
(or equivalently the Ri's) have their right-hand- sides that only depend on the data H and 
its partial derivatives, and on the Si 's to zero-th order. Then the work is to pass from this 
resulting system of PDE's to full gradient equations (i.e. total covariant derivatives) for 
the unknown basis {Si, . . . , Sn) or {Ri, . . . , Rn), equivalently. 

4.1 Elimination of source term and proof of formula ( fTHj ) 
4.1.1 The case n = 2 

First notice that because R G 5*02 (M), we have the relations 

JRi = e.jRj, {i,j) G h := {(1, 2), (2, 1)}, (34) 

where £12 = —€21 = 1 and we have defined J = ■ particular, this implies that 

V ■ i?, = [V, JRi] = e,, [V, Rj] , {i, j) G I2, 

as well as the relations, for any vector field A, 

[A, Ri] = JA- Ri = -A - JRi = -EijA ■ Rj, {i, j) G h- 

Together with the system of equations (!3T]) - (!32|) . these relations allow us to get the com- 
ponents of F in the basis {Ri,R2) (here, {i,j) G 12). 

F ■ Ri = —V ■ Ri + Vii ■ Ri + Vij ■ Rj 

= ~£^ij[V, Rj] + Vii ■ Ri + Vij ■ Rj 

= -e,j{[F, R,] + [Vji, R,] + [V,j, Rj]) + V, ■ Ri + Vj ■ Rj 

= —F ■ Ri + Vji ■ Rj — Vjj ■ Ri + Vii ■ Ri + Vij ■ Rj, 
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and hence the formula 

2F-Ri = -{Vii + Vjj) ■ + 2Vn ■ Ri + {Vi, + V^i) ■ Rj, (35) 

for {i,j) G l2- We then obtain equation f lTSj) by plugging f p5|l into the relation F = 
{F ■ Ri)Ri + {F ■ R2)R2, and using the following identity 

-V^ii-V^22 = Vlogrf, (36) 

whose proof may be found in section 14.1.31 below. 

4.1.2 The case n = 3 

We recall the vector calculus identity 'V-{AxB) = B- VxA — A-VxB, which holds 
for any two smooth vector fields A and B, together with the relations 

Ri{x) X R^{x) = R,{x), {i,j,k) G As := {(1, 2, 3), (2, 3, 1), (3, 1, 2)}, (37) 

which hold for all x G since R{x) G 5*03 (M). Combining these relations, we have 

V ■ R^ = V ■ {Rj X Rk) = Rk - V X Rj - Rj - V x Rk, {i, j, k) e A3. 

Using equations ( 1311) and ( l33l) . we obtain 

Vim ■ Rm — F ■ Ri = Rk ■ {Vjm X Rm + F X Rj) — Rj ■ (Vkm X Rm + F X R/.) 

= 2F ■ Ri — Vji ■ Rj + Vjj ■ Ri — Vki ■ Rk + Vkk ■ Ri- 

Above, we sum over m but not over (i, j, k) & A3. Thus we get 

3F ■ R, = -{Vu + V,, + Vkk) ■ + 2Vu ■ + {V, + V,,) ■ R, + {Vk + 14.) ■ Rk- 

We finally arrive at (1181) by plugging the three components F ■ Ri into the relation F = 
{F ■ Ri)Ri + {F ■ R2)R2 + {F ■ R3)R3 and using the following identity, 

_Vn-V22-V33 = V\ogD. (38) 

whose proof may be found in section 14.1.31 below. 

4.1.3 Proof of identities (|36l) and (|38l) 

We prove these identities in general dimension n > 2. First notice that for any k = 1 . . .n, 
the function T{x) satisfies locally the (tautological) relation 

dkT{x) = V\x)T{x), 

where we have defined V^ := {dk(tii)t'-^}ij = {dkT)T^^. Therefore by Liouville's formula, 
we obtain that 

n 

9fc(logdetT) = trace (V^'^) = ^1 

i=l 

Finally noting that detT = when n = 2 and detT = when n = 3, the identities 
( l36l) and ( 138|) are proved component by component. 
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4.2 Resolution in 5'0„(M) and stability analyses 

Plugging formula ( ITSj) into (13T1) . (l32l) and (!33l) provides a closed system of equations for 
the function R. It remains to show that R is then uniquely and stably determined by 
such equations. We again separate the cases n = 2 and n = 3. 

4.2.1 The case n = 2 

Since 5*02 (M) is a one-dimensional manifold, R = [i?i|i?2] is described by a §^-valued 
function 6{x) such that, at each point, Ri{6) = [cos 6 , sin 6)'^ and R2{d) = JRi{6). We 
wish to derive an equation for V6. Plugging the expression f fTSj) of F into fj3T|) . we arrive 

at 

V ■ i?. = ^ [-(Vlogd) ■ R, + {V,, - V,,) ■ R,] , (z, j) G /2. (39) 

Let us now derive a differential equation for R = [i?i|i?2]. The relation R^R = I2 
implies that R^diR G ^2(M), the space of two-dimensional anti-symmetric matrices. For 
z = 1, 2, i.e., it can be written in the form R^diR = aiJ, where ai = R2 ■ diRi. Defining 
OL := OiSi, it is clear that ct ■ Z = R2 ■ [{Z ■ \/)Ri] for any vector field Z. For the sequel, 
we need the following vector calculus identity, which holds for any smooth vector field A 

V|A|2 = 2(A- V)A-2[V,A]JA (40) 

We decompose a in the basis (-Ri,i?2) and use equations ( 1311) 

a = {a- Ri)Ri + {a ■ i?2)i?2 = {R2 ■ [{Ri ■ V)i?i])i?i + {R2 ■ [{R2 ■ V)i?i])i?2 

= (i?2 ■ [{Rl ■ V)i?l])i?l - (Rl ■ [(i?2 ■ V)i?2])i?2. 

We now use ( l40|) . which, in this case, becomes {Ri ■ V)Ri = [V, Ri]JRi = — (V ■ Rj)Rj for 
(i, j) G I2, and then equation (l39l) to obtain that 

a = -(V-i?2)i?l + (V-i?l)i?2 

= -^(-(1^12 - V21) ■ Rl - (Vlogd) • R2)Rl + ^(-(Vlogd) ■ i?i + (VS2 - V21) ■ R2)R2 

= ^[Vl2 - V21 - {R2 ®Rl-Rl® i?2)Vl0grf] = ^[Vi2 - V21 - JVlogrf]. 

Finally, given the above parameterization R{6), one checks using the chain rule that 
R^diR = {di9)J, i = 1,2. Using the preceding calculations, we obtain that 

V9 = a = ^ [V12 - V21 - JV log d]. (41) 

Reconstruction procedure. Let if be a given data set corresponding to an illumina- 
tion g that guarantees condition ( fT5|) throughout X. From 5", we construct R using the 
Gram-Schmidt procedure, and parameterize R with a function 6 : X ^ E>^ as above. Let 
Xm G dX be the global minimum of the illumination gi . Then according to the maximum 
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principle [TSj, Mi achieves its minimum over X at Xm- In particular, at this point, we have 
the relation 

{Xm) = -iy{x.m) = I'^^l^"!, = Rl{6{x„^)), 



|VMi| \Sl{Xjn)\ 

where ^{xm) denotes the outgoing normal vector to X at Xm- Thus 6{xm) is known. 
Therefore, we can reconstruct 6{x) at every x G X by integrating (HT]) along the segment 
] parameterized by 7x„,a;- 



■ / (^12 - V21 - JVlogd)(7,„,,(t)) dt. 

^0 



Once 6 is recovered throughout X, one then reconstructs cr{x) for all x G X from the 
knowledge of cr(xo) for some Xq G X and integrating equation (fT8|) along the segment 
[xo, x], that is 

loga(x) = loga-(xo) + 2 / j^oAt) " ^(7xo,x(t)) c^t 



= loga(xo) + (x - Xo) ■ / (Vlogrf + {{Vij + V,,) ■ i?,)i?,)(7,„,,(t)) dt. 

Jo 

Proof of Theorem 13.21 We now prove the stability results stated in Theorem 13.21 
Let two data sets {H,H) G W^'°°{X) correspond to identical illuminations g = g, and 
conductivities a and a that coincide at some Xq G X. Let T, T be built from H,H such 
that they satisfy conditions (125!) and (129|1 . We then define R = ST'^ and R = S'T^ and 
parameterize R and R by their angle functions ^, 6* : X — )■ as above. Using estimates 
and taking the difference of fHTl) for 6 and 6 yields 



2\V{e-9)\ < |V^2i-V^2i| + |^i2-V^i2| + |V(logd-logd)| <C\\H-H\\ (42) 

Moreover, if x^ G dX is the global minimum of the illumination gi, we have seen in the 
reconstruction procedure that R{6{xm)) = R{0{xm)) = —^-'{xm) and thus 6{xm) = 0{xm)- 
Therefore we have, for any x G X, 

^(x) - d{x) = d{x^) - ~d{xm) + {x-x^)- [ V{9 - ^)(7x„,.(t)) dt, 



and thus \9{x) - 9{x)\ < A{X)\\V{6 - ^)||l-(X), where A(X) denotes the diameter of 
X. Combined with fH2l) . this yields 

11^ — 6'||vFi'°°(x) < Cg\\H — vi/i,»=(x). (43) 

This inequality is also a stability statement for the reconstruction of the function 6. On 
to the stability of a, we have the pointwise estimate 

|V(loga - log a) I < |V(logrf - logrf)| + \{{Vpk - Vpk + Vkp - Vkp) ■ Rk)Rp\ 

+ \ {Vpk + Vkp) ■ {Rk — Rk)Rp\ + l((^pfc + Vkp) ■ Rk){Rp — Rp)\ 

< |V(logrf-logrf)| + 2 J2 \ypk-Vpk\ + M0-^\ l^p'^-l' 

l<k,p<2 l<k,p<2 
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where we have used the Cauchy-Schwarz inequahty and the fact that \dQRp\ < 1 for 
p = 1,2. Using estimates (143|) and (l29l) imphes directly 



||V(log(T - log(T)||L°o(X) < C'lll-f^ - H\\wi.<^{X)- 



The proposition is proved provided that a and a coincide at xo, which allows us to control 
any difference | \oga{x) — logcr(x)| by A(X)||V(log(T — loga)\\Loo(^x)- This concludes the 
proof. 

4.2.2 The case n = 3 

In the three-dimensional case, the unknown is now an 503(11^)- valued function R and is 
thus parameterized by 3 scalar functions in principle. However, because every 3-parameter 
chart on S'03(R) has singularities, as shown in [23], using such a chart to represent the 
unknown function R may yield instabilities during the reconstruction procedure. In this 
paper, we thus represent i? as a 9- vector R = , i?^, Rj)^. We then derive gradient 
equations for the components of R, which allow local reconstruction of R by integration 
of an appropriate system of ordinary differential equations. 

Remark 4.3. For numerical purposes, we want to choose a parameterization of SO^iM.) 
that involves as few parameters as possible while still guaranteeing stability. This can 
be achieved by the A-parameter quaternion representation, whose details will appear else- 
where. 

Let us now find differential equations for R. Taking partial derivatives of the rela- 
tion R^R = I3 yields that for k = 1,2,3, we have that R^dkR G ^3(M), the space of 
antisymmetric 3x3 matrices, which we write 



We now focus on expressing the vector fields ck, := a^^k in terms of the data. In order to 
do so, we recall that A3 is the set of direct permutations of (1, 2, 3). We notice that 



We now use the following identities, which hold for any smooth vector fields A and B 

2{B ■ V)A = V{A ■ B) + V X {Ax B) - Ax {V X B) - B X {V X A) - {V ■ B)A + (V • A)B, 
2{A ■ V)A = V|Ap -2Ax{V X A), 

where A = R^, B = Rp, and we take the dot product with Rj. For p = i, we obtain 

2cii ■ Ri = 2Rj ■ [{R, ■ V)Rk] 

= Rj ■ [V X Rj - X (V X R,) -Rix{Vx Rk) - (V ■ Ri)Rk + (V ■ Rk)Ri] 
= Rj ■ (V X Rj) -Ri-{Vx R,) +Rk-{Vx Rk). (46) 



R^ dkR 




k = 1,2,3. 



(44) 



ai- Rp = Rj ■[{Rp-V)Rk], {i,i,k)eA-i, pe{i,j,k}. 



(45) 
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For p = j, we get 



2cxi ■ Rj = 2R, ■ [{R, ■ V)Rk] = -2Rk ■ [{Rj ■ V)i?,] = 2Rk ■ [R^ x (V x 

= -2Ri ■ (V X Rj). (47) 



Finally for p = k, we get 

2ai ■ Rk = 2Rj ■ [{Rk ■ V)Rk] 



-2R.J ■ [Rk X (V X Rk)] = -2R, ■ (V x R^) 



(45 



We have just proved that the vector fields oli can be expressed in terms of the quantities 
{Rp ■ (V X Rq)}i<p,q<z- These quantities, in turn, can be expressed in terms of R and the 
data using equation (133|1 . For (i, j, k) G A3, we have 

(V X Ri) ■ Ri = {Va xRi + FxRi)-Ri = -V,j ■ Rk + Vik ■ R, 

(V X Ri) ■ Rj = {Vii xRi + FxRi)- Rj = Vu ■ Rk -Vik■R^ + F■Rk (49) 

(V X R,) ■ Rk = {Va xRi + FxRi)-Rk = -Va ■ R, + V, ■ R^ - F ■ Rj. 

To summarize, using (|T8|) for F and identity fl38|) . we have for k) G A3 that 
1 



ai-R^ = - [{Vkj - Vjk) ■ Ri - iVk + Vki) ■ Rj + {Vj, + Vij) ■ Rk] , 
■ = ^ + Vki) ■ R^ + {Vkj - 2Vjk) ■ Rj + (2V,,- + Vkk - Vii 



Rk 



(50) 



«i ■Rk = ^ [-{Vij + Vji) ■ Ri - {Vjj + 2Vkk - Vii) ■ Rj + {2Vkj - Vjk) ■ Rk] 
These calculations show that there exist vector fields Apqr such that 



CXr 



i^Apqy • Rq)Rj., 



I.e. 



Am Dm r)k 



k 



1,2,3. 



Left-multiplying equation (H4|) by i?, we obtain the differential equations 



(51) 





dkRi 




' -a^R2 + a^Rs ' 






dkR2 




-a\Rz + a^Ri 


, 1 < A; < 3 




dkRs 




-Q!2-Ri + "1-R2 





(52) 



Since the expression (I5T!) is purely quadratic in the components of R, the above system 
can be summarized as 



|a|=3 



(53) 



1=1 



where the functions : f2 — t- are linear combinations of the data vectors Vij ( l28ll . 

Local reconstructibility and stability estimates. The system (153|) allows us to 
reconstruct R locally provided that we know its value R(xo) = Ro for some Xq G Vt. The 
reconstruction of R(a:) for x G ^2 is done by integrating the following ODE 



d 
di 



R(7(t)) = G(7(t), R(7W)), R(7(0)) = Ro 



(54) 



16 



where 7 is chosen to be ^xo,x here, and where the function G is polynomial in the compo- 
nents of R with bounded coefficients provided that H has components in W^'°°. Because 
G is constructed such that G ■ R = 0, as can be seen in (!52|) . the solution to this system 
of equations, if it exists, has constant ||R|P norm, equal to ||Ro|P = 3. Therefore, the 
function R(a;) G remains in the compact set -\/3S^, over which G is Lipschitz in the R 
variable. Thus by virtue of standard results for ordinary differential equations [TB], the 
solution to fl5^ exists, is unique, and its existence can be extended up to t = 1. Such a 
reconstruction procedure is stable in the sense of the following proposition. 

Proposition 4.4 (Local stability when n = 3). Let Vt C X and H,H : Vt ^ S^^^R) he 
two data sets with components in W^'°°{VL) that satisfy the condition 

inf(D,5) > Co > 0. (55) 

Then there exist two constants Gq,Gi such that for every x,y we have the estimate 
||R - R\\{y) < CoWR - R||(x) + C^H - i/||iyi-(n). (56) 



Proof of proposition 4-4 Let two data sets H,H : Q ^ 5*3(11^) with components in 
W^'°°{Q) correspond to conductivities (a, a). Let T, T be built from H,H such that 
they satisfy conditions fl2S]) and fl2^ . We then define R = ST^ and R = ST^. For the 
rest of the proof, we denote by R = {Rj, R2 , Rj)^ and R similarly. Recall that we have 
equations of the type 

|a|=3 

where the Q^'s are linear combinations of {V^j} ( 128|) . Given the stability condition (129|) . 
we obtain estimates of the form 

max ||QQ||L°°m) < C'||-f^||vi^i.°°(n) and max \\Q^ — < C\\H — H\\wi,oc>m\. 

a,k a,k 

(57) 

We then write the pointwise relation 

dk\\R - Rf = 2(R - R) ■ dkiR - R) = 2(R - R) ■ ^ [Q^R" - g^R"] 

|a|=3 

|a|=3 

Recall that R and R satisfy ||R||^ = ||R||^ = 3 on so each of their components is 

lal = 3, 



bounded by Vs. Therefore, for any 9-index a with 

3 3 



|R" -R° 



1=1 1=1 



— K^U -Rji )-^i2 -^«3 I ~l~ Lli2)Ri^\ ~\~ \RiiRi2{Ri3 Ri. 

< 9||R- Rll. 
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Therefore we have the estimate 



a^llR - Rf < 2||R - R|| ^ [||Q^ - Q^ll^ + 9||gJ 

|a|=3 



|R- R| 



Writing c}fc||R — R|p = 2||R — R||9fc||R — R|| and using the estimates in (1571) . we obtain 
for some constants Co and Ci that 

(5fc||R — R|| < Co||i/ — i?||oo + C'i||iJ||oo||R — R||- 

Let x,y E Q. Then applying GronwalFs lemma for the function ||R — R|| over the segment 
[x,y] parameterized by ■-yx,y, we get the estimate 



||R - K\\{y) < (||R - R||(x) + AiQ)Co\\H - i7i|i,oo)e^^^(^)ll^"^-. 
This concludes the proof. 



(5^ 



□ 



4.2.3 Gram-Schmidt decomposition 

The Gram-Schmidt decomposition is defined as follows 



Ri := Si/\Si\ and Rj = [ Sj - ^(5^ ■ Ri)R, / 



1=1 



i-i 

Sj — ^(Sj ■ Ri)Ri 

i=l 



J>1, 



and builds a matrix R from S that is orthogonal by construction. In the three-dimensional 
case, the transition matrix T such that R = ST'^ is given by: 



{tij} 









_{H,2H23-H22H,s)idDy^ {H,2H,,-HuH23){dD)-^ dD~^ 



(59) 



with d := {HnH22 - and D = (det H)^ . 

The vector fields Vij defined in ( l28l) take the following expression 



1<«J<3 



V log til 

t22 ^7 *21 



til t22 

tu tz'i tiit22 tgg t22 h3 





Vlogt22 

V log ^33 



(60) 



When n = 2, T and Vij are given by the top-left 2x2 blocs of ( 1591) and ( 1601) . respectively 



Proof of lemma \4.1\ For simplicity we only prove the Lemma in the case where the di- 
mension is two. The three dimensions can be done similarly. We first notice that 
||-f^~^(2;)||2 < ll-f^ ""^lloo over Q. Since if is a symmetric matrix we have 



H-\x)\\2' < H,,{x), 
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and so \\H ^H^^^ < Hjj{x) over Q for 1 < j < 2 . Using these estimates and Remark 13.11 
we easily derive the inequahties 

\\tij — tjjlloo < IviIq'^ Mq^H — if 1 1 14/1,00 (Q), 
4 

i=o 

for 1 < i, J < 2, which proves the result. □ 



5 Global reconstructions in 3D 

As we have seen in the past section, the three-dimensional approach is locally similar to 
the two-dimensional far as the pointwise derivation of reconstruction formulas 

is concerned. The main difference is that, because the result [1] Theorem 4] does not 
hold when n = 3, one must work with a covering of X and build a global reconstruction 
procedure that patches the local reconstructions over the domains together in a stable 
manner. 



5.1 Proof of lemma 13.31 

Let us first justify that the approach presented in the next section is valid in the sense 
that we can build illuminations such that a couple {0,j) satisfies condition (l22l) . 

Proof of lemma \3.3[ A way to fulfill condition (l22|l is to construct solutions whose gradi- 
ents are taylored to satisfy this condition up to terms that can be made negligible. We 
do this by using the Complex Geometrical Optic (CGO) solutions, a generalization of 
the harmonic complex plane waves of the form e'''^ with p G C"' such that p ■ p = 0, so 
that Ae'''^ = p ■ pe^'^ = 0. These functions were first introduced in [11] in the context 
of linearized inverse problems, then extended in [25] in the context of non-linear inverse 
problems. The construction can be made for any n > 2. 

Construction of CGO's: We first extend the diffusion equation V ■ {a{x)'Vu) = to 
M", where (t{x) is extended in a continuous manner outside of X and such that a = 1 
outside of a large ball. Assuming that a\x G H2'^^~^'^(^X) for some e > (this implies 
a e C^{X) by Sobolev imbedding), it is shown in [9l Corollary 3.2] following works in 
[111 [25] . that there exist complex- valued solutions of the above full-space diffusion equa- 
tion of the form 

u, = ^e''-^{l + ^P,), (61) 



where p G C" is a complex frequency satisfying p ■ p = 0, which is equivalent to taking 
p = p(k + ik-'") for some k, k"*" G S"~^ such that k-k"*" = and p = \p\l \pl > 0. Moreover, 
the remainder ipp satisfies the PDE 

AV'p + 2p- V^p = — ^(1 + V^p). (62) 

'(J 
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Energy estimates on (1^ done in [H] use the regularity assumed on a to show that p'i)p\x = 
0(1) in C^{X\ Using this estimate, computing the gradient of (16T]) and rearranging term, 
we arrive at 

v/ffVup = e'''^(p + (^p), with ipp:=V^p ^ ^pp - {\ ^p)V ^fa . 

Because V^ct is bounded and p'^p\x = 0(1) in C^{X), the C"-valued function cpp satisfies 
supjf < C independent of p. Moreover, the constant C is in fact independent of a 
provided that the norm of the latter is bounded by a uniform constant in H2~^^~^^[X). 
Completion of the proof: Note that Up is complex-valued, thus its real and imaginary 
parts Mp and provide two solutions of the diffusion equation, and i/aVw^ and y/aVu^ 
may serve as vectors S^. More precisely, we have 

V^Vn^ = pe^^-^ ((k + V^) cos(pk^ ■ x) - (k^ + p-^?) sin(pk^ ■ x)) , 
v^Vm^ = pe"'^-" ((k^ + p-V?) cos(pk^ ■ x) + (k + p^Vp) sin(pk^ ■ x)) . 

Let us now define pi = p(e2 + iei), p2 = p{es + iei), and construct 
(5i, ^2, ^3, 54) = V^(V<, Vul, V<, Vul). 

Then, using (I63p . we obtain that 

det(5i, S2, Ss) = p3e''(2-2+x3) (_ + ^ 

det(5i, 52, 54) = p3e^(2x.+x3) ^_ ^-^^^^^^ ^ j^^^^^ ^ 

where limp^oo sup^ |/i| = limp_^oo supx I/2I = 0. Letting p so large that supx(|/i|, I/2I) < 
I and denoting 70 := min^gj^(p^e^*^^^2+2;3)-j ^ have that | det(S'i, 5*2, S's)! > ^ on 

sets of the form X fl {pxi g]^, |[+m7r} and | det(S'i, 5*2, 5'4)| > ^ on sets of the form 
X n {pxi g]|, ^[+m7r}, where m is a signed integer. Since the previous sets are open 
and a finite number of them covers X (because X is bounded and p is finite), we therefore 
have fulfilled the desired construction. Upon changing the sign of 5*3 or 5*4 on each of 
these sets if necessary, we can assume that the determinants are all positive. 

Let {gi]i<i<i be the traces of the above CGO solutions (m^^, m^^, m^^, m^J. These illu- 
minations generate solutions that satisfy the desired properties. Any boundary conditions 
Qi in an open set sufficiently close to gi will ensure that the maximum of the determinants 
stay bounded from below by Cq > 0: 

max(det(S'i, 6*2, 6*3), det(5'i, 6*2, S4))){x) > Cq > 0, x & X. 

This concludes the proof of the lemma. □ 



5.2 Global reconstruction procedure 

Let us consider m > 3 solutions of ( ITOl) and assume that there exist an open covering O = 
{f2j}i<j<Ar with X C UfLiQi, a constant Cq > 0, and a function r such that assumption 
(!22|) holds. Let us fix xq G fij^ such that cj(xo) and {Sr{iQ)i{^o)}i<i<3 known. We now 
assume that there exists an integer K > and two functions 

Y : X 3 x^ Y{x) = {yi{x) = xq, y2{x), . . . , yxix), yK+i{x) = x) e 

-.X x[l,K]3 {x,i)^il:{x,z) e[l,N], ^ ' 
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such that for every x G X, 

[xo,x] = Uf^^[yi{x),yi+i{x)] and [y^yi+i] C l<i<K. (65) 

An example of such a setting is shown in Figure [TJ 




Figure 1: An example of setting for the reconstruction of a at some x & X. One can 
see that [xo,y2{x)] C Qi, [y2{x) , y3{x)] C Q2 and [?/3(x),x] C 1^4. In this case, we have 
^{x, 1) = 1, i{x, 2) = 2 and ij{x, 3) = 4. 

From the 4-uple {0,T,Y,ip), we define the reconstruction procedure V{0,t,Y,iIj) as 
follows. Let X e X. For every 1 < k < N, \et S^^'^ = [5'^(^(x,fc))i |5'r(,/,{x,fc))2 |5'r(^(x,fc))3], 
f/'C') = {ifjjjj jg^(^(3, ^j)), and T'-'^^ the matrix of Gram-Schmidt coefficients of H^''\ For 
further use, we also define 

^k-i,k _ {Hij{yk)}i(ZT-{^(^x,k-i)),jeT{ij{x,k))- (66) 

The reconstruction procedure V{0,T,Y,tlj) is summarized in algorithm 15.11 
Here and below, the superscript ~^ denotes the matrix inverse transpose. That the 
transfer equation fl^7|l holds can be seen from the following facts: 

• one can compute S^'^^^^yk) from R^'^^^^yk) by writing 

S^'^-'\y,) = R^'^-'\y,) . {T^'-'\y,))-^. 

• Then, from the equality {S^^~^\yk)r S'^''\yk) = H''^^'^ with H^"^''' defined in (166]), 
and the fact that S^''~^\yk) is known, one obtains S^'^^yk) by writing S^'^^y^) = 

• Finally, one computes R^''\yk{x)) by formula R^^\yk{x)) = S^''\yk{x))T^''\yk{x)) . 

Combining the three formulas above gives (l67j) . 

In order for this procedure to yield globally stable reconstructions, we not only need 
local stability inside of each Qi ^ O (guaranteed by proposition 14.41) . but also stability 
as one passes from ^^(x,k-i) to ^^[x,k) using formula fl67|) . This is done in the following 
lemma. 
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Algorithm 5.1 Reconstruction procedure V{0,T,Y,tlj) 



for X G X and 1 < k < K do 

construct T^'^^ and the vector fields V^|^^^ (1251) from the measurements {Him}^m&T{ji>{x,k)), 
and define R^''^ = S'^'^-'T^^^-^ by orthonormalizing S^'^K 

if k = 1 then 

(xo) is known from the knowledge of S^'^\xo). 

else 

construct the value of R^^\yk) from the value R^^~^\yk) using the following formula 
R^^\y,{x)) = R^^-'\y^{x)) ■ T^^-'\yu{x)) ■ H'-''' ■ T^^\yu{x)). (67) 

end if 

For all y G [yfc(a;), compute R^^\y) by integrating flM|) along [yk{x),y]. 

Compute a{yk+i{x)) from (j{yk{x)) and R^''\y),y G [yfc(x), ?/fc+i(x)] by integrating 
equation ([18]) along [yk{x),yk+i{x)]. 
end for 



Lemma 5.1. There exist constants C2 and C3 such that for every x G X and every 
2 < k < K , we have 

||R(^)(y,(x)) - R^'\yk{xm < C2\\R^'-'\ydx)) - K^'-'\y,{xm + Cs\\H - H\U,^. 

(68) 

Proof. The transfer equation ( 1671) can be summarized as R^^\yk) = R^''~^\yk)M''~^''', 
where M''~^'^ := T^^~^\yk{x)) ■ H^-^''^ ■ T^^\yk{x)), and similarly for R^^\yk). Provided 
condition ([22]) holds, T^^'^^ and T^'^) are well-defined and we have estimates of the form 

llM'^-i'^ll < C||^||i,oo and ||M'^^i'^-M^-^'^|| <C||iJ-^||i,oo, (69) 
where ||M|| denotes any matrix norm. The lemma is thus proved if we write 

R^''\yk) - R'-'Kyk) = (R'-'-'Kyk) - R^^-'\yu))M''-^'^ + R!^^~^\yk){M''-'^^ - M'-^'"), 

and bound \\R'-''\yk) - R^^Kyk)\\ 

using triangle inequalities, estimates (l69l) and the fact 
that p('=-i)(yfc)|| < 1. □ 

We now conclude the proof of theorem I3.4[ 



Proof of theorem \3.4\ Let us define the error function 

e(x) := logcr(x) — log5"(x). 
Using equality (ITSl) . we obtain that 

^Ve{x) = V(logD(^) -log5(^)) 
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where i,j run over the set rliplx, K)) and the superscript (K) denotes the fact that we 
are at the i^-th (i.e. last) step of the reconstruction. By writing 

+ {{V,,j - V,,) ■ (R, - R,))Rj + {{Vij - V,,) ■ Ri){R, - Rj), 
we obtain the following bounds 

^||V£(a:)|| < C\\H - iJ||i,oo + 2||if ||i,oo||R(^H^) - R^''^^)!!- (70) 

We now apply alternatively K times proposition 14.41 and lemma 15.11 Define C4 = 
max(Co,C2) and C5 = max(Ci,C3) where Cq.Ci are given in proposition 14.41 and 02,0^ 
are given in lemma 15.11 We obtain 

||R(^)(x) - R(^)(^)|| < C4||R(^)(2/i^(^)) - K^^'KyKixm + C4H - E\\,^^ 

< Cl\\K^^-'\yK{x)) - R(^-i)(yA'(x))|| + C,{1 + C,)\\H - H\\,^^ 

< ... 

<C|^||R«(yi(x))-R«(yi(x))|| 

+ C,{l + C,■■■ + Cl''-')\\H-H\\^^^ 

_ f<2K _ 1 _ 

= Cn|R«(xo) - R^'\xo)\\ + Cs-^ -\\H - H\U,^. (71) 

04 — i 

Using the fact that R*^^^(xo) and R'^-'^)(xo) are built from {S'^(jg);(xo), 5'T-(jQ)j(xo)}i<i<3, we 
obtain an inequality of the form 

3 

||R(i)(xo) -R«(xo)|| < CrJ2\\Sj.J^o)-S,^Jxo)\\+Cs\\H-H\\i,oo. (72) 

1=1 

Summing up the three bounds ( !70l) . ( ITTi) and (172|) . we obtain that 

3 

||Ve(x)|| <CJ2\\Sri^o)X^o) - Sri^o),{Xo)\\ + C'\\H - H\\,,^. 
1=1 

Considering that e{x) is bounded by |e(a;o)| + A(X) || V£:||oo for any x E X, inequality ( j24l) 
is proved. □ 
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